#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <math.h>
#include <stdbool.h>
#include <float.h>

#include <quadmath.h>
#include <immintrin.h>

static inline __float128 i32_to_float128(int i)
{
   char numbuf[64];
   sprintf(numbuf,"%d",i);
   return strtoflt128(numbuf, NULL);
}

void max_exp_table_prec52()
{
   __float128 q1 = 1.0Q;
   __float128 q2 = 2.0Q;
   __float128 q52 = 52.0Q;
   __float128 q128 = 128.0Q;

   __float128 r1;
   __float128 x;
   __float128 y;

   printf("#include <stdint.h>\n\n");
   printf("// exp_data_52[i]= (2^(i/128)-1.0)*2^52\n");
   printf("uint64_t exp_data_52[128]=\n{\n");
   for (int i=0;i<128;i++)
   {
      r1 = i32_to_float128(i);
      x  = r1 / q128;
      y  = (powq(q2,x)-q1)*powq(q2,q52); // y=  (2^(i/128)-1.0)*2^52
      y  = roundq(y);
      printf("  0x%013" PRIx64 ,(uint64_t)(y));

      if ( i!=128-1)
         printf(",");
      if ( i % 2==1)
         printf("\n");
   }
   printf("};\n");
}

void max_exp_table_prec64()
{
   __float128 q1 = 1.0Q;
   __float128 q2 = 2.0Q;
   __float128 q64 = 64.0Q;
   __float128 q128 = 128.0Q;

   __float128 r1;
   __float128 x;
   __float128 y;

   printf("#include <stdint.h>\n\n");
   printf("// exp_data_64[i]= (2^(i/128)-1.0)*2^64\n");
   printf("uint64_t exp_data_64[128]=\n{\n");
   for (int i=0;i<128;i++)
   {
      r1= i32_to_float128(i);
      x = r1 / q128;
      y = (powq(q2,x)-q1)*powq(q2,q64); // y=  (2^(i/128)-1.0)*2^64
      y = roundq(y);
      printf("  0x%016" PRIx64 ,(uint64_t)(y));

      if ( i!=128-1)
         printf(",");
      if ( i % 2==1)
         printf("\n");
   }
   printf("};\n");
}

void max_exp_table_prec61()
{
   __float128 q2 = 2.0Q;
   __float128 q61 = 61.0Q;
   __float128 q128 = 128.0Q;

   __float128 r1;
   __float128 x;
   __float128 y;

   printf("#include <stdint.h>\n\n");
   printf("// exp_data_61[i]= 2^(i/128)*2^61\n");
   printf("uint64_t exp_data_61[128]=\n{\n");
   for (int i=0;i<128;i++)
   {
      r1= i32_to_float128(i);
      x = r1 / q128 + q61;
      y = powq(q2,x);
      y = roundq(y);
      printf("  0x%016" PRIx64 ,(uint64_t)y);
      if ( i!=128-1)
         printf(",");
      if ( i % 2==1)
         printf("\n");
   }
   printf("};\n");
}


int main(int argc, char* argv[])
{
#if defined(PREC_52)
   max_exp_table_prec52();
#elif defined(PREC_61)
   max_exp_table_prec61();
#elif defined(PREC_64)
   max_exp_table_prec64();
#else
   #error must define PREC_52 or PREC_61 or PREC_64
#endif
  return 0;
}
